#### "Economic Inequality. Income and political participation in authoritarian regimes" ####
# authors: "Pelke, Lars"
# date: 2019-06-17
# written under "R version 3.6.0 (2019-03-11)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(haven)
library(lme4)
library(sjstats)
library(texreg)
library(interplot)
library(ggrepel)
library(stargazer)
library(ggpubr)
library(dotwhisker)
library(broom.mixed)

#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))

#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A068>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 |A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 | A102>=0
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 )%>%
  mutate(civilsociety_b = if_else(A064==1, 1,
                       if_else(A066==1, 1,
                       if_else(A067==1, 1, 
                       if_else(A068==1, 1,
                       if_else(A069==1, 1, 
                       if_else(A070==1, 1,
                       if_else(A072==1, 1, 
                       if_else(A073==1, 1,
                       if_else(A074==1, 1,
                       if_else(A075==1, 1,
                       if_else(A076==1, 1, 
                       if_else(A077==1, 1,
                       if_else(A078==1, 1,
                       if_else(A079==1, 1, 
                       if_else(A099>=1, 1,
                       if_else(A100>=1, 1, 
                       if_else(A101>=1, 1,
                       if_else(A102>=1, 1,
                       if_else(A103>=1, 1,
                       if_else(A104>=1, 1, 
                       if_else(A105>=1, 1,
                       if_else(A106>=1, 1,
                       if_else(A106B>=1, 1, 
                       if_else(A106C>=1, 1, 0)))))))))))))))))))))))))

table(wvsdata$civilsociety_b)
mean(wvsdata$civilsociety_b)

number_observations <- wvsdata %>%
  count(country_year)

# Dependent variable without party organizsations 
wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 | A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 ) %>%
  mutate(civilsociety_p = if_else(A064==1, 1,
                          if_else(A066==1, 1,
                          if_else(A067==1, 1, 
                          if_else(A069==1, 1, 
                          if_else(A070==1, 1,
                          if_else(A072==1, 1, 
                          if_else(A073==1, 1,
                          if_else(A074==1, 1,
                          if_else(A075==1, 1,
                          if_else(A076==1, 1, 
                          if_else(A077==1, 1,
                          if_else(A078==1, 1,
                          if_else(A079==1, 1, 
                          if_else(A099>=1, 1,
                          if_else(A100>=1, 1, 
                          if_else(A101>=1, 1,
                          if_else(A103>=1, 1,
                          if_else(A104>=1, 1, 
                          if_else(A105>=1, 1,
                          if_else(A106>=1, 1,
                          if_else(A106B>=1, 1, 
                          if_else(A106C>=1, 1, 0)))))))))))))))))))))))

table(wvsdata$civilsociety_p)

#Replacing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$civilsociety_b)

# 40217 non member of civil society/ not belonged to civil society, 32846 member/ Belong to cso

number_observations2 <- wvsdata %>%
  count(country_year)

## Inspect which countries are delete because missing obserations at macro-level variables ##

setdiff(number_observations$country_year, number_observations2$country_year)

# Comment: [1] "China 1990"        "Croatia 1996"      "Georgia 1996"      "Hong Kong 2014"    "Jordan 2007"       "Mexico 1981"      
# [7] "Mexico 1990"       "Russia 1990"       "South Africa 1982" "South Korea 1982" are excluded from the dataset

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 

### delete Kuwait 2014, Uzbekistan 2011 and Libya 2014, no GINI information in a range of maximum 4 years before or after the country-year observations

wvsdata <- wvsdata %>%
  drop_na(gini_mkt)
rm(missing)

table(wvsdata$civilsociety_b)
# 36955 non member of civil society/ not belonged to civil society, 30677 member/ Belong to cso

#### Previous democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 


### delete missing macro-level cases 

wvsdata <- wvsdata %>%
  drop_na(e_migdppcln)

#### Country and Country-Year overview ####

country_year_observations <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_year_observations
rm(country_year_observations)

country_observations <- wvsdata %>% 
  group_by(country_name) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_observations
rm(country_observations)

count <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  tally() 
rm(count)

number_observations3 <- wvsdata %>%
  count(country_year)
setdiff(number_observations2$country_year, number_observations3$country_year)

## Deleted observations because of missing macro-level variables: 
# [1] "Kuwait 2014"         "Lebanon 2013"        "Libya 2014"          "Palestine/Gaza 2013" "Uzbekistan 2011"    
# [6] "Zanzibar 2001" 

rm(number_observations, number_observations2, number_observations3)

##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene), 
         regime_support_group_max = ifelse(is.na(regime_support_group_max), 0, regime_support_group_max))

# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2csreprssCGM <- centFUN(wvsdata$v2csreprss)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)

summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2csreprssCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)


#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
wvsdata$ageCWC <- wvsdata$ageCWC/10

##################################
#Descriptive Staticstics
##################################

country_observations <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_observations
rm(country_observations)

count <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  tally() 

wvsdata <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  mutate(civil_byincome = mean(civilsociety_b))

wvsdata <- merge.data.frame(wvsdata, count)

rm(count)

ggplot(data=wvsdata, aes(x=income_quintiles, y = civil_byincome/n))+
  geom_col() +
  facet_wrap(~country_year, ncol = 5) +
  theme_bw() +
  xlab("Income Quintiles") +
  ylab("Civil Society Participation") +
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Average Participation in civil society by income groups",
       x = "Income Quintiles",
       title = "",
       caption = "Data from World Value Survey and V-Dem dataset, version 10")

ggsave("output/A2.pdf", width = 22 , height = 25 , units = c("cm")) # Figure B2 in Supplementary Appendix

#### Scatterplot Inequality and Civil Society Participation ####

library(ggrepel)

wvsdata_country <- wvsdata %>%
  group_by(country_year) %>%
  summarize(voting_civil = mean(civilsociety_b),
            inequality = mean(gini_mkt))

ggplot(data= wvsdata_country, aes(x=inequality, y = voting_civil, label = country_year)) +
  geom_point() +
  geom_smooth(method = lm, colour = "black") + 
  geom_text_repel(size = 3) +
  theme_classic() +
  labs(y = "Average Civil Society Participation in Country-Year",
       x = "Market Inequality",
       title = "",
       caption = "Data from World Value Survey and V-Dem dataset, version 10")

ggsave("output/A1_scatter_civil_inequality.pdf")

rm(wvsdata_country)

#### Histogramm Civil Society Participation ####

ggplot(data=wvsdata, aes(x = civilsociety_b)) +
  geom_histogram(binwidth = 1, fill = "white", colour = "black") +
  theme_classic() +
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("No", "Yes")) +
  labs(y = "Number of individuals",
       x = "Civil Society Participation (yes/no)",
       title = "",
       caption = "")
ggsave("output/A1_histogramm_civil.pdf")

##Summary statistics ##

cols <- c("sex", "ageCWC", "educationCWC", "children", "unemployed", "income_quintilesCWC",
          "gini_mktCGM", "gini_dispCGM", "e_migdppclnCGM", "v2csreprssCGM", 
          "regime_support_group_maxCGM", "democratic_expericeneCGM", "civilsociety_b" )

library(stargazer)
stargazer(
  title="Summary Statistics for Civil Society Particiaption Dataset", 
  covariate.labels=c("Male","Age (centered)", "Education (centered)", "Children", "Unemployed",
                     "Income Quintiles","Market Gini (c)","Disp. Gini (c)" ,"CSO repression (c)",
                     "GDP pc log (c)", "Regime Support Group", "Democratic Experience", 
                     "Civil Society Participation"),
  as.data.frame(wvsdata[, cols]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd")
)


##################################
#Analysis: three-level models: Market Inequality
##################################

v1 <- glmer(civilsociety_b ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
icc(v1, adjusted = TRUE)
icc(v1)


v2 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)

#macrolevel gini
v3 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)
isSingular(v3, tol = 1e-4)


#macrolevel gini plus addtional macro controls
v4 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)
isSingular(v4, tol = 1e-4)


# model 4 plus random slopes
v5 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM +   
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)


# model 5 plus interactions 
v6 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*gini_mktCGM + v2csreprssCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)
isSingular(v6, tol = 1e-4)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2csreprssCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)
isSingular(v7, tol = 1e-4)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_mkt <- interplot(v6, 
                                       var1 = "gini_mktCGM", 
                                       var2 = "income_quintilesCWC", 
                                       ci = 0.95, 
                                       hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_civil_mkt
ggsave("output/C1interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(v6, 
                                          var1 = "v2csreprssCGM", 
                                          var2 = "income_quintilesCWC", 
                                          ci = 0.95, 
                                          hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C1interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(v7, 
                                          var1 = "v2csreprssCGM", 
                                          var2 = "income_quintilesCWC", 
                                          ci = 0.95, 
                                          hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C1interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors",  "Market Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_main_civil <- {dwplot(v_models,
                           dot_args = list(size = 2, aes(shape = model)),
                           whisker_args = list(size = 0.5, aes(linetype = model)),
                           dodge_size = 1) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.75),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                    name = "Model", 
                    breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                    labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                         labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"))
  
  } %>% 
  add_brackets(v_brackets)

plot_main_civil
ggsave("output/C1dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))

plot_main_vote <- {dwplot(v_models,
                          dot_args = list(size = 2, aes(shape = model)),
                          whisker_args = list(size = 0.5, aes(linetype = model)),
                          dodge_size = 1) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                      labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                         labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"))
  
} %>% 
  add_brackets(v_brackets)

plot_main_vote

#### texreg regression tables main models ####

library(texreg)
texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:C3", 
       longtable = TRUE)

###############################
##Using Disposable Inequality##
###############################

## Supplementary Appendix: C2 ##


d1 <- glmer(civilsociety_b ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)
icc(d1, adjusted = TRUE)
icc(d1)


d2 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)
isSingular(d4, tol = 1e-4)

# model 4 plus random slopes
d5 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus interactions 
d6 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*gini_dispCGM + v2csreprssCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)


# model 6 but with interaction of freedom/fariness and income level and without gini
d7 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2csreprssCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d6, d7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_disp <- interplot(d6, 
                                        var1 = "gini_dispCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_disp
ggsave("output/C2_interaction_inequality.pdf")


plot_Interaction_civil_repress <- interplot(d6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_quintilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C2_interaction_repressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(d7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_quintilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C2_interaction_repressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

d_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

d_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Disp. Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_main_civil_disp <- {dwplot(d_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(d_brackets)

plot_main_civil_disp
ggsave("output/C2_dot_whisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:C4", 
       longtable =  TRUE)